      !<@file
      !>@brief Subrotina para a entrada dos parâmetros de inversão.
      !>
      !> Todas os parâmetros necessários para solução do problema inverso são inseridos por meio desta rotina: algoritmo de otimização, tipo de rgularização, parâmetros dos algoritmos, critérios de parada, modelo inicial, modelo de referência, etc. 
      subroutine entradainversao(maxIter,step,tol,limit,eps,xtol,msave,cmin0,cmax0)

      use globais
      implicit none

      double precision:: cminref,cmaxref !< valores extremos para criação do modelo de referência com variação linear vertical
      double precision:: cmin0,cmax0!< valores extremos para criação do modelo inicial com variação linear vertical
      double precision :: STEP  !--STEP SIZE ALONG SEARCH DIRECTION FOR FINAL ITERATION, parâmetro do CG   
      double precision :: Tol  !--COMPUTING TOLERANCE (ITERATIONS STOP WHEN MAX-NORM OF GRADIENT .LE. Tol) , parâmetro do CG    
      integer :: MaxIter! parâmetro do CG (não é o número máximo de iterações)   
      integer :: limit!< número máximo de iterações, parâmetro do CG   e do LBFGS
      integer :: Hnnzmax!< número máximo de elementos não nulos da matriz H esparsa(aproximação da Hessiana)
      integer :: i,j,kk,ndonoh
      integer :: it !< contador de instantes de tempo

      integer::ierr

      double precision :: eps !< determina a acurácia da solução, parâmetro do L-BFGS
      double precision :: xtol !< tolerância da busca linear, parâmetro do L-BFGS
      integer :: msave !< número de correções usada no update do BFGS, parâmetro do L-BFGS

      real, dimension(:,:), allocatable :: Pauxreal !< variavel auxiliar para leitura dos sismogramas
      double precision, dimension(:,:), allocatable :: Pauxdouble !< variavel auxiliar que recebe os valores lidos como real, convertendo para double (motivo: A rotina dft precisa de double!)
      complex*16, dimension(:), allocatable :: Pwaux !< variavel auxiliar para a transformada dos sismogramas
      character*5 :: chks !< caracter para o número da fonte
      character*2 :: chigsr!< caracter para o número do grupo de pares fonte/recptor

      allocate(PreceivBARRA(Nrec(1),Nsource(1),Nomega,gruposSR),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif












      open(1,file='parametros.txt',status='old')

      
      read(1,*)Nx0, Nz0
      N=Nx0 * Nz0
      read(1,*)cminref,cmaxref
      read(1,*)cmin0,cmax0

      read(1,*)metodo,nloops,clust,suav
      if(metodo.eq.1)then !< CG
	  read(1,*) STEP, tol !=0. ! STEP =0:  o algoritmo escolhe o passo inicial da busca na direcao!= 1.0d-6
	  !MaxIter = n ! sugestao do algoritmo (num max de iteracoes ateh renormalizar a direcao)
	  MaxIter = 0 ! renormaliza a cada iteracao (jah que eu faco poucas iteracoes, muito menos q o numero de incognitas)
	  !write(*,*)'maxiter=', maxiter
      elseif(metodo.eq.2)then
	  read(1,*) eps, xtol, msave
      endif
      read(1,*) limit! = Numero máximo de iterações (nome proveniente do CG)

      read(1,*) Kreg, lambdanaux ! tipo de regularizacao
!       lambdanaux = lambdan
      close(1)

      if(kreg.eq.4 )then
	  qreg = lambdanaux
      endif

      deltaX0 = Lx0/dfloat(nx0)
      deltaZ0 = Lz0/dfloat(nz0)

      write(*,*)'# Razao entre delta do PI e do PD (dir x e z):',deltax0/delta,deltaz0/delta

      write(*,*)'alocando C0...'
      allocate(c0(nx0,nz0),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando Xref...'
      allocate(xref(n),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando CXr...'
      allocate(CXr(n),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif




      write(*,*)'alocando resid...'

      allocate(resid(Nrec(1)*Nsource(1)),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif




      write(*,*)'alocando residexpandido...'
      allocate(residexpandido(Nx*Nz,Nsource(1)),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando gradmisfit...'
      allocate(Gradmisfit(N),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando treg...'
      allocate(GradTreg(N),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando gradFobjetivo...'
      allocate(GradFobjetivo(N),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando deltaR...'
      allocate(deltaR(Nx*Nz,Nrec(1)),stat=ierr) ! du/dc(x)  vou usar esta matriz fonte por fonte      
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando JtdeX...'
      allocate(JdeX(Nrec(1),N),stat=ierr)     ! du/dX
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif
      

!       Hnnzmax=N ! se Hdiagonal
!       Hnnzmax=5*N ! se H tiver 5 diagonais, hessiana de hu
      Hnnzmax=N*N ! se H completa

      write(*,*)'alocando Hvalue...'    
      allocate(Hvalue(Hnnzmax),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando Hrowind...'    
      allocate(Hrowind(Hnnzmax),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando Hcolind...'    
      allocate(Hcolind(Hnnzmax),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando Hcolptr...'    
      allocate(Hcolptr(N+1),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      write(*,*)'alocando Hprec...'    
      allocate(Hprec(N,N),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif

      if(cminref.lt.cmin  .or.  cminref.gt.cmax)then
	  write(*,*)'erro no modelo de referencia'
	  stop
      else if(cmaxref.lt.cmin  .or.  cmaxref.gt.cmax)then
	  write(*,*)'erro no modelo de referencia'
	  stop
      else if(cmin0.lt.cmin  .or.  cmin0.gt.cmax)then
	  write(*,*)'erro no modelo inicial'
	  stop
      else if(cmax0.lt.cmin  .or.  cmax0.gt.cmax)then
	  write(*,*)'erro no modelo inicial'
	  stop
      endif

      do j=1,nz0
	  do i=1,nx0
	      ndonoh = (j-1)*nx0+i
	      Xref(ndonoh) = cminref  + (cmaxref-cminref)/dfloat(nz0-1)*(j-1)
	  enddo
      enddo

      do i=1,N
	  Xref(i) =  (  Xref(i) - (cmax+cmin)/2  )*2.2/(cmax-cmin)!isso eh uma aproximacao da funcao nao linear de X em C por uma reta
      enddo


      ! armazena quais nohs do grid cada incognita influencia
      sizedep = (int(deltaX0/delta)+1+Npml)*(int(deltaZ0/delta)+1+Npml)


      write(*,*)'alocando depX...'
      allocate(depX(sizedep,N),stat=ierr)
      if(ierr.eq.0)then
	  write(*,*)'.........Ok'
      else
	  write(*,*)'.........allocation error',ierr     
	  stop
      endif



      do j = 1, N
	do i=1,sizedep
	    depX(i,j) = 0
	enddo
      enddo
      call criadepCX()




      if(kreg.ne.0)then
	  ! plota modelo de referencia se houver regularizacao
	  write(*,*) 'Modelo de referencia'
	  call modelo2d(Xref)

! 	  call systempsimage(c,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'veloR',&
! 	    & 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
! 	    & wbox,Lz0/4.d0,Lx0/4.d0,cmin+1000.d0,cmax-1000.d0,(cmax-cmin-2000.d0)/4.0d0, &
! 	    & 0.0d0,Lz0,0.0d0,Lx0,1)!modelo do pd sem PML

! 	  call systempsimage(c,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'veloR',&
! 	    & 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
! 	    & wbox,Lz0/4.d0,Lx0/4.d0,cmin+1000.d0,cmax-1000.d0,(cmax-cmin-2000.d0)/4.0d0, &
! 	    & 0.0d0,Lz0+(Npml-1)*delta,-(Npml-1)*delta,Lx0+(Npml-1)*delta,1)!modelo do pd com PML


! 	  call systempsimage(c0,nx0,nz0,deltaX0,deltaZ0,deltaZ0/2.0d0,deltaX0/2.0d0,'veloR',&
! 	    & 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
! 	    & wbox,Lz0/4.d0,Lx0/4.d0,cmin+1000.d0,cmax-1000.d0,(cmax-cmin-2000.d0)/4.0d0,&
! 	    & 0.0d0,Lz0+(Npml-1)*delta,-(Npml-1)*delta,Lx0+(Npml-1)*delta,1)!modelo do pi, sem a pml

      elseif(kreg.eq.0)then
	  call system('rm veloR.bin veloR.ps &')
      endif


      ! Entrada dos dados


! 	! allocate(norm2PBARRA(Nomega))
! 	open(1,file='medidasfreq.txt',status='old')
! 	read(1,*)!Titulo   '#Nreceptor Nfonte Nfrequencia Medida'
! 	do kw=1,Nomega
! 	    do igsr=1,gruposSR
! 		do ks=1,Nsource(igsr)
! 		    do kr=1,Nrec(igsr)
! 			read(1,*)PreceivBARRA(kr,ks,kw,igsr)
! 		    enddo
! 		enddo
! 	    enddo
! 	enddo
! 	do igsr=1,gruposSR
! 	    do ks=1,Nsource(igsr)
! 		do kr=1,Nrec(igsr)
! 		    write(*,*)(PreceivBARRA(kr,ks,kw,igsr),kw=1,Nomega)
! 		enddo
! 	    enddo
! 	enddo
! 	close(1)

      call system('tar -xzf sismD.tgz') 
      write(*,*)'sismogramas descompactados!'



      write(*,*)'Carregando dados dos sismogramas...'
      do igsr=1,gruposSR ! Loop no número de grupos fonte/receptor

	    allocate(Pauxreal(nttotal,nrec(igsr)),stat=ierr)
	    allocate(Pauxdouble(nttotal,nrec(igsr)),stat=ierr)
	    allocate(Pwaux(nomega),stat=ierr)
       

	  do ks=1,Nsource(igsr) ! Loop nas fontes

	      call loadSismog(Pauxreal,nttotal,nrec(igsr),igsr,ks,'sismD')


	      do kr=1,nrec(igsr)
		  do it=1,nttotal
		      Pauxdouble(it,kr) = DBLE(Pauxreal(it,kr))! isso pq a rotina dft foi feita para double precision
		  enddo
	      enddo


	      do kr=1,Nrec(igsr)
		  call dft(Pwaux,omegas,Nomega,ntps,Pauxdouble(1,kr),nttotal)
		  do kw=1,Nomega
		      PreceivBARRA(kr,ks,kw,igsr)=Pwaux(kw)
		  enddo
! 		    write(*,*)(PreceivBARRA(kr,ks,kw,igsr),kw=1,Nomega)
	      enddo

	  enddo !< loop das fontes



	  deallocate(Pwaux)
	  deallocate(Pauxreal)
	  deallocate(Pauxdouble)

      enddo!< loop dos grupos de pares fonte/receptor




      return

      end subroutine entradainversao